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ABSTRACT 



We consider vertical heat transport in Keplerian accretion disks, including the 
effects of radiation, convection, and turbulent mixing driven by the Balbus-Hawley 
instability, in astronomical systems ranging from dwarf novae (DNe), and soft X-ray 
transients (SXTs), to active galactic nuclei (AGN). In order to account for the interplay 
between convective and turbulent energy transport in a shearing environment, 
we propose a modified, anisotropic form of mixing-length theory, which includes 
radiative and turbulent damping. We also include turbulent heat transport, which 
acts everywhere within disks, regardless of whether or not they are stably stratified, 
and can move entropy in either direction. We have generated a series of vertical 
structure models and thermal equilibrium curves using the scaling law for the viscosity 
parameter a suggested by the exponential decay of the X-ray luminosity in SXTs. We 
have also included equilibrium curves for DNe using an a which is constant down to a 
small magnetic Reynolds number (~ 10^). Our models indicate that weak convection 
is usually eliminated by turbulent radial mixing, even when mixing length estimates 
of convective vertical heat transport are much larger than turbulent heat transport. 
The substitution of turbulent heat transport for convection is more important on 
the unstable branches of thermal equilibrium S-curves when a is larger. The low 
temperature turnover points ^max on the equilibrium S-curves are significantly reduced 
by turbulent mixing in DNe and SXT disks. However, in AGN disks the standard 
mixing-length theory for convection is still a useful approximation when we use the 
scaling law for a, since these disks are very thin at the relevant radii. In accordance 
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with previous work, we find that constant a models give almost vertical S-curves in 
the S — r plane and consequently imply very slow, possibly oscillating, cooling waves. 



Subject headings: accretion, accretion disks - novae, cataclysmic variables - X-rays: 
stars 



1. Introduction 

The disk-instability model in a thin accretion disk is the most popular model for dwarf-novae 



(DNe) outbursts and soft X-ray transients (SXTs) (for recent reviews see Cannizzo 1993 or Osaki 



19961) . In this model the initial fast rise of luminosity in an outburst is due to the propagation of a 
heating front at the local thermal speed, while the exponential decay follows from the subsequent 
slow propagation of a cooling wave. The alternation of heating and cooling in these systems 
comes from the characteristic 'S' shaped curve in the midplane temperature versus column density 
equilibrium curve (see Fig. |l| below). At all points along this curve the heating due to the 
local dissipation of orbital energy is balanced by radiative losses and convective heat transport. 
The unstable middle branch of the equilibrium curve, where column density falls with increasing 
midplane temperature, is produced by the hydrogen ionization transition and its effect on radiative 
opacities. The local heating rate is usually taken to follow some variant of the a prescription 



( Shakura and Sunyaev 19731 ). Vertical heat transport in accretion disks is modeled following 



the usual treatment for radial heat transport in stars, with radiation and convection defining 
the thermal energy balance. Unlike main-sequence stars which have stable convective zones and 
internal veloicty fields, convection in accretion disks varies with time as they evolve along the 
thermal equilibrium curves. Typically, convection has a very strong influence on the structure of 
the thermally and viscously unstable middle branch of the thermal equilibrium curve since the 
optically depth of a partially ionized disk is very large. 

Convective mixing is not the only source of turbulent transport in stars and disks. Diffusive 
mixing induced by shear instabilities in rotating stars has been examined by several researchers 
(Zahn 1992, Maeder 1995, Maeder and Meynet 199^ ) for its effects on the luminosity, and evolution 



of stars. Similarly, it is appropriate to include shear- induced turbulent mixing in the study of 
accretion disk structure, especially since the luminosity of accretion disks is ultimately due to 
the dissipation of orbital energy through some sort of shear instability. Turbulent mixing is an 
unavoidable part of this process. Furthermore, a particular local instability has been identified 
(Velikhov 1959, Chandrasekhar 1960, Balbus and Hawley 1991, Hawley and Balbus 1992) which 
is uniquely suited to playing this role, and consequently gives us a theoretical basis for evaluating 
the efficiency of turbulent mixing. In accretion disks the Velikhov-Chandrasekhar instability 
is commonly referred to as the Balbus-Hawley instability. Here we abbreviate this as the BH 
instability. In contrast to convective transport, which only affects the convectively unstable regions 
and moves heat down temperature gradients, turbulent transport driven by the BH instability is 
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ubiquitous within magnetized accretion disks, and will always act to reduce entropy gradients. 
The effect of turbulent diffusion due to some sort of shearing instability, not necessarily the BH 
instability, has been studied in the context of radial transport in thin (important for transition 
fronts, e.g. see Menou, Hameury, & Stehle 1999) and advective ( |Honma 1996D disks, and in the 



context of vertical transport in T Tauri disks ( D'Alessio et. al. 199^ ) 



In a real disk the BH instability should generate turbulence which is only moderately 



anisotropic ( Vishniac and Diamond 1992 ), an expectation which is consistent with the available 
numerical simulations (see Gammie 1997 and references contained therein). Consequently, we 
expect that the dimensionless form of all the transport coefficients should differ only by factors of 
order unity. More precisely, we note that the numerical simulations consistently suggest that the 
angular momentum transport is dominated by magnetic tension, and that the relevant magnetic 
field moment < B^Bq > is a few times p < >. The vertical diffusion coefficient is < > t, 
where r is the velocity correlation time scale, and Jlr is a factor of a few, where is the local 
orbital frequency. We conclude that the vertical diffusion coefficient is comparable to the radial 
diffusion coefficient. In other words, 

Dz ~ V, (1) 

where v is the usual notation for the turbulent kinematic viscosity. Here we will assume 
that equation (Q) is an exact equality, although in reality we expect that this is only a rough 
approximation. Consequently, if we assume that the turbulent Prandtl number is unity, the 
vertical thermal diffusion coefficient xt can be written as 

where Q is the Keplerian angular velocity, and P and p are the disk pressure and mass density. 
The turbulent mixing rate across the disk is 

klxT ~ a^^, (3) 

where we have assumed a vertical wavelength l/k^ comparable to the vertical global pressure scale 
height Ap. Since this is also the thermal equilibration rate for the disk, this implies that turbulent 
mixing can have a significant effect on disk structure. 

The BH instability follows from the presence of a magnetic field in the disk, which requires 
some sort of magnetic dynamo. The value of a is a fraction of (Va/c^)^, the ratio of magnetic 
pressure to ambient pressure. This ratio is determined by the balance between turbulent dissipation 
and flux generation via the (uncertain) dynamo process. It has long been known that the observed 
luminosities and timescales for DNe require a larger hot state viscosity, ahot^ than in quiescence, 
o-coid (Smak 1984, [Livio and Spruit 1991 ). This could indicate a dependence of a on temperature. 



or more probably on the disk geometry. On the other hand, it could simply reflect the fact that 
cold disks have a large neutral fraction and a high resistivity. The latter consideration has led 
Gammie and Menou (1998) to propose that disks in quiescence are cold enough to shut down all 
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dynamo activity. Evidence in favor of the former idea comes from the observed exponential decay 
of soft X-ray luminosity for soft X-ray transients which seems to require a scaling law for a of the 
form ( [Cannizzo, Chen, and Livio 1995 , Vishniac and Wheeler 199(: . Vishniac 1997) 



««ao(^) , (4) 

where ao = 50, Cg is the isothermal sound speed at the mid-plane, and r is the radial position of 
the annulus. The value of ao is proportional to the assumed masses of the central black holes 
used to calibrate this relation. Cannizzo et al. assumed a black hole mass of IOMq. This scaling 



law has been criticized (see e.g. King &: Ritter 1998 and references therein) on the grounds that 



irradiation should substantially modify the structure and stability of disks around black holes. 
Certainly a naive calculation of the optical emission from the disks of SXTs during the decline 
from outburst fails to match the observed behavior ( |Chen, Shrader, and Livio 1997| ). However, 
a recent study ( Dubus et al. 1999| ) has shown that it is much more difficult for self-consistent 



models of irradiation to change the basic disk structure (and by extension the X-ray light curve) 
than to alter the optical emission from the disk. 

Equation (^) is not compatible with the original suggestion by Balbus and Hawley (1991) 
that the dynamo process should saturate in equipartition between the magnetic field and the gas 
pressure. This suggestion was motivated by the notion that the dynamo process arises from the 
random stretching of field lines in a turbulent medium ( Batchelor 1950| , Kazantsev 1967). On the 



other hand, this kind of dynamo is rigorously justifiable only in the limit where magnetic forces are 
negligible, which is obviously not the case here. Equation (^) may be compatible with the notion 
that the instability gives rise to an incoherent dynamo effect ( [Vishniac and Brandenburg 1997 ), 



although this predicts a somewhat steeper dependence on geometry as (cs/ril) 0. However, 
the saturation of the instability as the magnetic field approaches equipartition may explain the 
discrepancy ( Vishniac 1999| ), and is, in any case, necessary to explain the low sensitivity to 



geometry seen in the numerical simulations (whose parameters are generally equivalent to Csjr of 
order unity). We stress that since the dependence is fundamentally one of geometry rather than 
Mach number, the current generation of numerical simulations are compatible with either model. 
Here we will explore the consequences of equation (Q) for vertical structure, although we will also 
discuss the alternative, that a is constant down to a limiting value of the magnetic Reynolds 
number. 

Another major issue in vertical structure calculations is the role of convection. Although 
convection is unlikely to account for any significant level of radial angular momentum transport 
(oc {vrVg), where Vr and are the fiuctuating velocity fields) in Keplerian disks ( [Balbus and 



Hawley 199^ ), it may still play an important role in vertical energy transport (oc (v^dT), where 



6T is the fluctuating temperature field). The use of standard mixing length theory (MLT) is even 
more uncertain in disks than in stars, and previous structure calculations have included simple 
nonlinear closure schemes (Canuto and Hartke 1986, pabot et. al. 198% pannizzo 1992 ), detailed 



linear analyses of the convective modes ( [Ruden, Papaloizou, and Lin 198^ ), and a reduced effective 
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length annizzo 19921) . Some of these calculations included the assumption that angular 



momentum transport is mediated by convective motions, so that the results cannot be extended to 
ionized disks. There is an additional problem with the use of nonlinear closure schemes since such 
schemes typically fail to reproduce experimental and numerical results in extremely strong shearing 
environments, although there has been some progress in the intervening years (cf. Speziale et al. 
1991; Speziale and Gatski 1997). A particular problem with such schemes is that they tend to 
underestimate the degree of anisotropy for low Rossby numbers. In any case, none of the previous 
work included the effect of diffusive damping and mixing due to turbulent motions acting across 
the reduced radial length scales. How can we determine the 'typical mode' of convection (i.e. 
in the spirit of MLT) when differential rotation and BH turbulence are simultaneously present? 
To account for nonlinear effects within the context of linear theory, we start with the linear 
perturbation theory of convection for Keplerian disks and find the modes most resistant to the 
dissipative effects of shear and turbulent mixing. We use the properties of these linear modes as a 
guide to the nonlinear state and the typical velocity and eddy size of convective cells. We defer 
discussion of the empirical basis for this approach to the last section of this paper. 

The purpose of the paper is to calculate the vertical structure and thermal equilibrium 
S-curves for a variety of astrophysical disks, including the effects of turbulent energy transport 
throughout the whole thickness of the disk, reasonable modifications of MLT for convection in 
disks, and the scaling formula for a given in equation (Q). Standard MLT and the associated 
vertical energy transport normally used in Keplerian disk models are reviewed in §2. Section 3 
describes the turbulent energy flux driven by the BH instability. This energy flux is particularly 
important if convection is suppressed. We examine the conditions for the suppression of convection 
in §4. In §5, we propose a new model for vertical energy transport in disks. We present and discuss 
our results for DNe, SXTs, and AGNs in §6. Finally, we summarize our results and present our 
conclusions in §7. 



2. Standard MLT with a scaling law for a 

The basic equations for the vertical structure of a geometrically thin, Keplerian, optically 
thick disk are derived from mass conservation, radiative and convective energy transport, and 
thermal and hydrostatic equilibrium. These give four differential equations for the column density 
E, the total vertical energy flux F ^ the total pressure P, and the temperature T as a function 



of distance from the disk midplane: ( Meyer and Meyer-Hofmeister 1982 , C annizzo and Cameron 



19881) 



f = = (7) 
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where z is the vertical coordinate perpendicular to the disk plane and is zero at the mid-plane, /x 
is the mean molecular weight, k is the Rosseland opacity, Frad is the radiative flux, V and are 
defined as d\nT/d\nP associated with the radiative flux F^ad and the total flux F respectively 



(|Cox and Giuh 196j 



4ac ^ 

Frad = pV. (9 

F = ^^V.. (10) 
The gas pressure Pg, is the difference of the total pressure and radiation pressure; i.e., 

Pg = P- \aT\ (11) 

We note that V^, the gradient d\nT/d\nP when the total energy transport is carried by radiation, 
does not necessarily correspond to any real gradient. It is a mathematical convenience used to 
calculate the total flux F. In our model the turbulent energy flux is always nonzero, although 
sometimes negative, and the total flux is never carried by radiation alone. We will address the 
nature of the turbulent energy transport later. 

These four differential equations have four boundary conditions. The deflnition of the column 
density T, = 2 J pdz gives S(z = 0) =0. The reflection symmetry about the mid-plane gives 
F{z = 0) = 0. The other two boundary conditions follow from the deflnition of the photosphere, 
P{z = h) = \gz{z = h)/K{z = h) and F{z = h) = aT^{z = h), where h is the half-thickness (from 
the mid-plane to the bottom of the photosphere) of the disk. 

We need one more differential equation for the variable a. Since we lack a model for the 
vertical distribution of dissipation, we will assume that the scaling law for a depends only on the 
mid-plane sound speed. This gives us a trivial differential equation 

da ^ , . 



with the boundary condition 



where 



a(z = 0) = 50('^^^^)''', (13) 



There are good reasons to think that equation ( |l^ ) is not correct (see, for example, flgure 10 in 
Stone et. al. 199"^ ), and will have to be modifled once we understand the z dependence of the 



dynamo process. 
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Now we turn to the coohng processes which play important roles in disks in the relevant 
temperature range, i.e. thousands to tens of thousands of degrees K. The usual approximations 
for vertical energy transport in a Keplerian, optically thick accretion disk include the Rosseland 
approximation for the radiative energy flux Fj-ad and MLT for the convective energy flux Fconv.MLT- 
Hence, the total vertical flux in the convectively unstable region is given by 

4ac \ 1 

Frad + Fconv,MLT = ^ + -pVMLT^CpT— (V - V) , (15) 

where vmlt is the vertical speed of the convective bubble in the mixing-length theory, A is the 
mixing length for convection, Cp is the specific heat at constant pressure and V is dlnT/dlnP 
estimated inside the convective bubble. Notice that using V' instead of the adiabatic gradient 
V(jd implies heat exchange between a convective bubble and its environment over the course of its 
lifetime. This tends to erase entropy inhomogeneities and reduce the convective flux. In general, 
this process is described by the equation 

r' V - V 

where F', the convective efficiency, is given by 

8acT3 ■ ^ ' 

Here we have assumed isotropic convection so that the volume to surface ratio is V/A = A/6. The 
convective speed in this expression is vmlt = — 

iKN /2y/2, the square of Brunt- Vaisala frequency 
with radiative damping iV^ = —gzQ{S ~ V')/-^p> and rj' is the ratio of the excess heat generated 
within the convective bubbles to the excess energy radiated by the convective bubbles during their 
lifetime, and Q is the coefficient of expansion. Here the energy exchange could be the radiative 
losses given by F', or it could include any other diffusive effect given by r/'. We shall see that this 
is one of the crucial points in this paper. Equation ([l5| ) can be rewritten in terms of Vr and the 
convective efficiency for convective bubbles V . 

= V + aor'(V - V), (18) 

where qq = 3. Equations (^) and ( [T8| ) were originally used for stellar models and their use for 
disks is somewhat questionable. The important feature of traditional MLT is implicit in our 
expression for V which describes nearly isotropic convection without any drag on the motion of 
convective bubbles[|. We shall see below that equation (|l6|) should be modified to take into account 
the unique properties of accretion disks: shear as a result of differential rotation and turbulent 
mixing, and turbulent drag, due to the BH instability. 



^To be more precise, traditional MLT assumes that the turbulent drag is due only to the interaction of convective 
bubbles with each other, and can be represented by an efficiency for the conversion of the free energy into the kinetic 
energy of the bubbles. 



- 8 - 



3. Turbulent Heat Flux Driven by the BH Instability 

As discussed above, it seems reasonable to assume that viscous lieating and angular 
momentum transport in hot accretion disks are driven by the BH instability. Just as convective 
mixing leads to vertical energy transport through mixing, turbulent eddies driven by BH instability 
will also produce vertical energy transport at rate proportional to the local entropy gradient. 
Hereafter the term "turbulent flux" will be used to denote the energy flux induced by the BH 
instability. As noted before, in principle the turbulent thermal conductivity is not necessarily 



equal to the turbulent viscosity ( Riidiger et. al. 1988 ). However, given the approximately isotropic 



nature of the BH instability and for the reasons given in the discussion of equations (||) and 
we will take xt = ^- By analogy with the formula for the convective flux, we can obtain the 
turbulent flux by replacing the convective diffusion coefficient vmlt-^ with z^. In other words, the 
turbulent flux driven by BH stability is given by 

Fturb = IpuCpT^iV - V"), (19) 

where V" is (dlnT/dlnP) estimated inside a turbulent bubble driven by the BH instability. 
Introducing V" means that we also assume that there is an exchange of radiative energy between 
turbulent eddies driven by BH instability and their surroundings. We can define a turbulent 
efficiency T" for turbulent eddies analogous to T' for isotropic convection. In other words, replacing 
the term (vmlt^) with in the expression for T' , we have 

~ Sac r3 ' 

and its relation with the entropy inhomogeneities associated with turbulent cells is analogous to 
that for convective cells: 

r" V - V" 
_^ = 1 (21) 

where r]" is the ratio of the excess heat generated within the turbulent bubbles to the excess 
heat radiated from the bubbles during their lifetime. A similar set of equations for the radiative 
damping of the turbulent energy flux due to shear instabilities in a differentially rotating star 
has been studied by Maeder (1995). In physical terms, the radiative efficiency is the ratio of the 
radiative cooling time scale to the eddy turnover time. The cooling time scale due to radiation in 
the optically thick limit is the square of the radiative diffusion length associated with the turbulent 
bubbles divided by the radiative conductivity AacT^ / {?>K,p^Cp). 

There are several reasons to believe that turbulent energy transport could be very important. 
First, convection will vanish when entropy per mass increases away from the midplane, while 
the turbulent flux driven by BH instability is always present. Moreover, the turbulent flux can 
carry energy toward hot mid-plane when a disk is convectively stable. Second, the turbulent flux 
inferred by equation (|l^) has its maximal value when V' is equal to V ad- In this particular case of 
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weak radiative losses, equation ( [Tg| ) is roughly equal to 

^(V-V,d)~^, (22) 

Ap As 

where As is the entropy scale height. The ratio of this flux to the radiative flux (~ Mil^) is 
~ Xp/Xs, which will be of order unity unless a disk is nearly adiabatic ( [Vishniac 1993| ). In other 
words, the turbulent energy flux can be as large as the total energy flux when radiative losses are 
small. Finally, and most importantly, convection can suffer from strong shearing in a Keplerian 
disk and from strong mixing in a turbulent disk. We expect that turbulent mixing will be 
important whenever convective mixing is suppressed in spite of a unstable entropy gradient. We 
address this point in the next section. 



4. The Suppression of Weak Convection 

Modeling convection is often considered one of the main sources of uncertainty in stellar 
models. Despite this MLT has provided a reasonable framework for achieving an approximate 
understanding of stellar interiors and has been broadly applied to accretion disks. However, the 
question of its applicability to accretion disks has to be raised since strong shearing and mixing 
can be expected to have a dramatic effect on the disk environment, and to vitiate the usual 
assumption of nearly isotropic convection. Shearing will shorten the radial width of convective 
cells whenever convection is weak (i.e. the Brunt- Vaisala frequency —iN « Q). Turbulent 
mixing may prevent the formation of the kind of organized convective cells assumed in MLT. Here 
we will show how to modify MLT to take these important effects into account. 



4.1. The radial deformation of convective cells 

We begin by considering the linear perturbation equations for nearly incompressible 
(Boussinesq) , axisymmetric convection in an accretion disk some local viscosity. They are 

(iu) + k'^u)vr + ikr^i/ - 2VLV0 = 0, (23) 
{iuj + k^v)ve + ^Vr = 0, (24) 

{iu) + k'^v)v:, + ik^-^ + 5gz = 0, (25) 
/pi/r\ 

iuj6 - v^d^ In f — — j = 0, (26) 

krVr + kzVz = 0, (27) 

where ui is the frequency measured by a local observer rotating with the disk, 5 = Sp/p is the local 
fractional density perturbation, ^ = 6P/p is the pressure perturbation divided by the density, Vr, 
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ve, and Vz are velocity perturbations, kz and kr are the vertical and radial wavenumbers, and the 
vertical gravitational field is given by Qz = zQ,'^ + 2ttG'E{z), where the second term is from the 
self-gravity of the disk. Self-gravity is not negligible for AGN disks for small a ( ^akimoto and 



Coroniti 1981, |Lin and Shields 198(: , Cannizzo 1992| ) which, as we shall see later, can happen when 



we use equation (^). This system of equations implicitly assumes that the vertical wavenumber 
is small compared to the pressure scale height, and other vertical scale lengths. Of course, this 
is wrong, but it makes only a small difference if we use these equations and identify kz with the 
inverse pressure scale height at the end. 

A more serious issue is that we have simply ignored the inevitable loss of axisymmetry in 
the convective turbulence. Again, this is not as bad as it first appears. The existence of strong 
shear will tend to suppress high angular wavenumber modes. Physically this follows from the 
fact that modes with short angular wavelengths will respond to strong shearing by deforming 
into a spiral wave with a short radial wavelength. This will have an overwhelming effect on 
mode structure when the shearing time is shorter than the growth time of an unstable mode. 
Complete suppression of the instability follows when the radial wavelength becomes longer than 
the characteristic scale for the variation of the comoving frequency. In other words, we require 

ke « kr^, (28) 

where 7 is the convective growth rate measured by a local observer rotating with the disk, and 
kg is the azimuthal wavenumber. When kg satisfies this constraint, the effects of shearing will be 
modest. In general, this will make kg no larger than kz- When we calculate the volume to surface 
ratio for a typical convective cell shape below we will assume that kg ^ kz- 

Finally, we have left the effects of turbulent damping out of the entropy equation (equation 
(p^ . We do this to avoid double counting, since we add this effect to the usual effect of radiative 
damping when calculating the efficiency of convection. 

We can get a rough understanding of how the local shear modifies the shape of convective 
cells by considering the case u = 0. Then equations ([2^) through (p^) imply 



2 



rg-7 



A, (29) 



where A = kl/k"^, Tc = —iN = —iQzQiy — ^ad)/^p, and 7 = —iCo is the convective growth rate. 
If convection is weak in the sense that Tc « ^, A « 1 and the radial length of convection is 
shortened by the strong shear. In fact, we have A ~ (Tc/ri)^ if the dominant mode has a growth 
rate which is some fraction of Tc- 

The radial thinning of convective cells under the influence of the local shear has important 
consequences for the suppression of convection. First, this radial distortion makes it much easier 
for turbulent diffusion, due to the BH instability, to damp out convective cells. The damping 
rate due to radial mixing becomes vkz/A ~ vA^'^/A which is stronger than the isotropic result, 
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V / ^ aO. (assuming that A ~ Xp). Using our full perturbation equations p3| ) - ( p7|) we obtain 
the dispersion relation 

/„ l2,,\2\ / /T \ 2 , 1,2, A /„ l2,,\ 

(30) 

We see that a sharp increase in k'^ /k"^ due to the effects of shear is associated with much smaller 
values of 7. 

The second effect associated with the eddy distortion is that the average radiative diffusion 
length is now reduced. This means that the expression of radiative efficiency equation (17) needs 
to be modified for radially deformed convective cells. The definition of the radiative efficiency 
gives (|Cox and Giuli 1968|) 

where It is the distance between the center of the convective element, where the temperature 
perturbation is a maximum, and its surface. The ratio V/A is the volume to surface ratio for a 
convective bubble. Under normal circumstances It is nearly proportional to V/A and is A/2 for 
the isotropic case. We can make a plausible estimate for this length scale by comparison with 
V/A = A/6 in the isotropic case. We use It = 3(y/A) for an anisotropic parallelepiped and 
assume for simplicity that the convective cells have a parallelepiped shape. Then the expression 
for radiative efficiency becomes 

^ ~ 8acT^ ^' ^'^^> 
where we have defined the anisotropy factor 

/^[3(y/^)/(A/2)f = ^^^^^. (33) 

This has a maximum value of one when A = {kz/kr)'^ = 1, that is, in the isotropic limit. Note that 
this reduction in convective efficiency has nothing to do with turbulent mixing, and is purely the 
result of shearing. 



Equation ( |3^ ) shows that the radiative efficiency is reduced when the eddies are slow, 
Vz < vmlt, due to the turbulent drag and radially thin, that is / << 1. We have already 
incorporated turbulent drag into the dispersion relation (cf. equation (pO[)). Radiative damping 
can be added by replacing the Brunt- Vaisala frequency for the adiabatic case with one that uses 
V', as discussed in section 2. The dimensionless gradients can be related to T' (equation (|32|) ) 
through equation (|T6|). Before we do this, we need to include one last damping effect, the role of 
turbulent mixing in smoothing entropy gradients between convective bubbles. This is the topic of 
the next section. 
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4.2. Mixing losses 

The diffusion, and eventual loss, of the entropy inhomogeneities associated with convective 
bubbles due to turbulent mixing is similar to the effects of radiative loss ( p?alon and Zahn 1997 



D'Alessio et. al. 1998D and can therefore be treated in an analogous fashion. We do this by 



extending equation (|16| ) which describes efficiency of convection in the face of radiation transfer 
between entropy inhomogeneities. The diffusive term rf in equation ( [T^ ) can have either sign, 
depending on the physical context. For example, convection in extremely hot stellar interiors can 
involve a local competition between local thermonuclear energy generation and neutrino losses 



(Cox and Giuli 1968|). In a similar way, entropy transfer through turbulent mixing in accretion 



discs can be modeled by defining rl as the ratio of the turbulent to radiative heat transfer between 
convective bubbles. The turbulent mixing loss under pressure equilibrium is equal to the turbulent 
conduction (l/2)pi/CpAT, where AT is the temperature gradient within the convective bubbles. 
The radiative loss is 4acr^AT/(3K/9) in the optically thick limit. The ratio between these two. 



using definition of F" given in equation (20), gives 

rl = -aoF", (34) 

which gives a decrease in convective efficiency. 

In summary, shearing and turbulence act together to reduce convective heat transport. 
Shear alone shortens the radiative diffusion length, turbulent mixing alone smooths away the 
entropy inhomogeneities associated with convective bubbles, and the interplay of these two gives 
a turbulent drag on convective eddies. We can account for these three factors by modifying 
equations (16) and ([T7|) using equation ([3^) to obtain 



Cp kp^VzA J, f 1 ^ V - V 



1-7?' 8ac r3 "^U + aoF"j V - V,/ ^^^^ 

When these corrections lead to a suppression of convection, the turbulent heat flux induced by BH 
instability can have a significant effect on the vertical structure of accretion disks. As we shall see 
below, weak convection, where the convective growth rate is less than $7, is typically suppressed in 
a wide variety of accretion disks. However, before we can show this we need to determine 'typical' 
convective eddy sizes and speeds allowing for shearing and turbulent mixing. We do so in the next 
section. 



5. A New Model for Vertical Energy Transport 

A full treatment of vertical energy transport needs to include radiation, convection, and 
turbulent mixing driven by BH instability. As discussed above, MLT needs to be corrected by 
incorporating the effects of shear and turbulent mixing. The starting point for our corrections is 
the linear dispersion relation equation ( |30|) for convection. However, it would be unreasonable to 
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where 



expect the use of purely linear perturbation equations to give us an accurate model of the fully 
nonlinear end state. The linear dispersion relation gives a family of solutions, any one of which 
is an equally valid starting point, but whose properties can vary tremendously. This means that 
we can give physically reasonable estimates for the properties of the nonlinear regime from the 
linear theory only if we are able to find a 'typical' linear mode whose properties are the best guide 
to the nonlinear state; i.e., the mode which is most resistant to viscous damping and secondary 
instabilities. We are not going to propose a model for the turbulent spectrum for convection. 
Instead, our objective is to produce estimates for the typical vertical velocity Vz and the typical 
radial extent of a convective cell (i.e. still one-mode model), that reduce to the usual mixing length 
estimates when convection is strong and small scale turbulence is weak, but which incorporate 
physics neglected by MLT. 

The growth of the convective instability is ultimately limited by dissipative turbulence 
generated from secondary instabilities and from the interactions between convective modes. We 
can quantify this by considering the ratio of the mode shear to the growth rate squared, or 

Ps = — , (36) 

T 

Sij = ^{vi,j -Vj^i), (37) 

and 7 is the mode growth rate. When this is large the mode will quickly succumb to micro 
instabilities, so by minimizing Pg we can choose the linear mode most relevant to the nonlinear 
limit. This criterion is not, by itself, adequate because we are also dealing with the effects of 
turbulent mixing on small scales, which will introduce an effective viscosity into the system. This 
suggests that we also need to minimize the ratio of viscous damping to the mode growth rate. 
Since this is an independent effect, we will add the two in quadrature by defining a second ratio 

P..^-!^ (38) 

T 

and minimize 

Pu + Ps- (39) 

The relative weighting of these two terms is not obvious, but an equal weighting seems the most 
reasonable choice. Note that when the system is heavily damped, will be large for all modes, 
and minimizing the sum will reduce to a minimization of P,^. This guarantees that this criterion 
will always pick the modes that can still grow as we approach the limit of marginal stability. 

For an axisymmetric, incompressible convective roll, we have 

4S' = k\l + ^vl (40) 
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We can combine this with the expression for vg from equation (^), replacing iio with 7, to obtain 

^.^fJ:±4e. + ^U^V. (41) 



16(7 + A;2i/)2 4^ y V 7 

The last factor of this equation should be approximately independent of mode shape (A). The 
dispersion relation equation ( pO[ ) in the inviscid, non-shearing, and isotropic limit gives 7 = Fc/ \/2. 
Since we need to reproduce MLT in this limit we have 

V, = V2vMLT^ = A7/2. (42) 

In MLT we take = 2/A in MLT, that is, the perturbation wavelength is half the size of an 
isotropic eddy. Consequently, the last factor in equation (^) is equal to one. In the absence of 
rotation and viscosity, minimizing Pg implies = k1 or = A. In this limit Pg does not depend 
on the viscosity, and the influence of viscosity on mode shape is felt entirely through Py. 

We note that our treatment assumes that the vertical scale length k""^ remains the same as 
in standard MLT, A/2. This follows from our use of the minimization principle. Both Ps and 
Pij involve wavenumbers through the ratio of k^ to k^, which does not affect the value of /c^, and 
through the strength of viscous effects, which always favor a small k^. Therefore we assume that 
k~^ is as small as possible, about one pressure scale height. 

With this in mind, we need to minimize 



tot 



16(7 + A;2zy)2 4^/ 7 



as a function of A and 7, with 



,2 8 l + AXpg, 

This minimization is subject to constraints from the dispersion relation equation (3C) and to the 
total energy transport equations involving radiation, convection, and turbulent mixing driven by 
BH instability. The purpose of the paper is to solve the differential equations (^), (^), (|^), (|8[), 
and (|l^) with the corresponding boundary conditions to investigate the vertical structures and, 
furthermore, to see how the disk structure changes when one uses the scaling law for a. 

To solve these equations, we need to express the variable V (or say Fj-ad) in the equation (^) in 
terms of the thermodynamical variables we want to solve. Since the formulae for V depend upon 
the vertical energy transport mechanisms, the thermal gradient V has different forms within and 
outside convection zones. Furthermore, the usual criterion for convective instability by comparing 
Vr with Vad needs to be modified since V,- is never physically relevant, due to the constant 
presence of turbulent flux. Therefore, we adopt a convective criterion involving the turbulent flux, 
where the thermal gradient V,- is replaced with Vrc which is (dlnT/dlnP) of the medium when 
convection is absent ( p'Alessio et. al. 1998 ). Using this new criterion, we are ready to solve for V 



inside and outside the convection zones. 
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We begin with the case rc > "^ad- Then all three energy transport mechanisms coexist 
within the convectively unstable region. Hence the total vertical energy flux, F, is 

F — Prad ~\~ Pconv ydamp ~\~ Pturb (^^) 

4oc 1 1 , „ 1 1 , 

Combining this equation with the definition of V,. gives 

V, = V + aor'i (V - V) + aor" (V - V") , (46) 



with the definition of V in equation ( |3^ ) and the definition of Vz is equation (^). Moreover, the 
convective efficiency V and the turbulent efficiency T" can be related to the Vs by using equations 



and (21). Since the heat generation by viscosity is proportional to aP, which is a constant 
inside and outside the convective bubbles, r]" should be zero and ??' = —aor" as discussed above. 

Combining equations (^), (|3^), (pSj), ( p2[ ) and ( |46| ) we can find the thermal gradients in 
terms of F', T", / and therefore in terms of 7 and A. We obtain 

V = (1 - C)V. + CV,d, (47) 

where 

1-C.fl+ +^)". (48) 



Note that Q has lost its original meaning referring to convective efficiency ( Cox and Giuli 1968| ) 
and used here purely for convenience. Equation (^7|) is one constraint to the minimization problem 
arising from the total energy transport. Using it we can write Tc as 



= ^(V-V) (49) 

9zQ r' 



(i-C)(v.-v 



ad) 



Xp F + l + aoF' 

Now we can combine all the constraints together by substituting this expression for into the 
dispersion relation, so that the result becomes the only constraint on the minimization problem in 
terms of the two variables 7 and A. Once the minimization is solved numerically, the value of C is 
determined by 7 and A. This in turn gives V in terms of and Vad using equation (^7|). We get 
the ratio of the radiative flux to the total flux {F^ad/F = V/Vr) in the differential equation (P) 
and then solve the four differential equations. 

Outside the convection zone, V rc < ^ ad, we have 

= V + aoF" (V - V") , (50) 

F" - (51) 

2 
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Solving for V, we have (Maeder and Meynet 1996) 



V = pTT fv^ + Vad] ■ (53) 



We can see the physical meaning of the turbulent efficiency in this case (Maeder and Meynet 



19961 ). When T" goes to infinity, we have 



V7 V7 V7" V7 -^turb V Vad 

v-v,,, V -V,,, = - V, • 

Note that when turbulent transport is very efficient, it does not necessarily imply that the 
turbulent flux carries of the flux. When T" goes to zero, we have 

V^V., V"-V., E^ = Yi_X^Yi_1l = o, (55) 

r V r Vr 



Finally, V in equation (|50|) gives the expression of Vrc ■ 

Vr + aor"V' .... 

Vrc is equal to V in equation (|5^) outside the convection zones. We note that using Vrc 
instead of does not matter much in the end. When they differ, it will be because of either 
strong turbulent mixing or weak convection. In either case turbulent mixing will tend to suppress 
convection altogether. 



6. Results and Discussion 

The computer program used to solve the vertical structure equations is a modiflcation of 
the one described in two previous papers ( pannizzo and Wheeler 1984 , Cannizzo and Cameron 



1988). The adiabatic gradient Vad, mean molecular weight /x, and coefficient of expansion 
Q = —din p/dlnT were supplied by H. Saio for these papers using the input physics given in 
Iben 1963| . We assumed a solar abundance of Na, K, Fe (assuming that the overall ionization 



potential xm = 5.76 eV), and C. We use the Rosseland opacities n of |Cox and Stewart 1969 for 



logT{K) > 3.9, those of [Alexander, Johnson, and Rypma 1983 for log T(i^) < 3.9, and those 



of pollack, McKay and Christofferson 1985| for T{K) < 1220. We smoothed logK linearly over 
3.8 < logT{K) < 4.0 and 1220 < T{K) < 1675 to ehminate the slight discontinuity in those 
temperature ranges as done in the original code. In the calculation, we set the mixing length A for 
convection equal to max(Ap, z). The nonlinear minimization problem for the subjective function 
Ptot subject to the equality constraint equation (|30|) with positive variables 7 and A is achieved 
by the Fortran routines D0NPL2 ( [Spellucci 199^ ) and NLPQL QSchittkowski 1985| ). 



Fig. m shows the thermal equilibrium S-curves for a DN at r = lO^'''^ cm with a white dwarf 
mass = 1 Mq. Every point in the thermal equilibrium curves is one solution to the differential 
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equations (||) to (|8|) and each curve is associated with a different model for the vertical energy 
transport. We do this to illustrate the effect of the different mechanisms and to show how our 
model differs from the previous standard model. Here we show the results of a model with 
heat carried by radiation and MLT convection {Frad + Fconv^MLx), by radiation and turbulent 
flux driven by BH instability {Frad + Fturb), by all three transport mechanisms but with MLT 
convection {Ftot,nodamp), or by all three transport mechanisms with modified (damped) MLT 
[FtotAamp)- We have also added a fifth curve showing the sensitivity of our results to the choice of 

"0- 

These curves differ mainly on the middle, unstable, branch because the disk is partially 
ionized at these temperatures and the resulting opacity peak lowers the efficiency of radiative 
transfer relative to convective or turbulent mixing. The surface densities at the left and right 
turnovers of the unstable branches are usually denoted by S^m and Smax- In this case, only S^ax 
shows any sensitivity to the input physics. The short stable branch within the unstable branch, 
which appears in some of the curves, is an effect of very strong mixing. It has been previously 
shown, for convective mixing, that the unstable branch moves up and to the right in the T — S 
plane as the mixing length A increases (Mineshige and Osaki 1983, pojmaiiski 1986 , Ludwig, 



Meyer-Hofmeister, and Ritter 1994 ). Increasing the mixing length is equivalent to enhancing 
the strength of convection. This suggests that regardless of the mixing mechanism at work, the 
upper-right shift and the intermediate branch of the unstable branch indicate the presence of 
stronger mixing. This explains why both the standard MLT convection model and the addition 
of turbulent mixing without damping produce qualitatively similar curves, with the latter lying 
consistently to the right. On the other hand, the curve for radiation and turbulent mixing alone 
lies on top of the curve with damped convection. This joint curve lies consistently to the left 
at lower temperatures, although the old standard model (F^ad + Fconv,MLT) crosses it near the 
intermediate branch. The fact that the addition of damped convection produces no significant 
changes indicates that turbulent damping is very efficient in this case. The tendency of the joint 
curve to lie below, and to the left, of the othe curves shows that vertical heat transport via 
turbulent mixing is generally less effective than undamped convective transport. The difference 
between the joint curve and the old standard model largely vanishes at high temperatures 
because the value of a, and the consequent strength of turbulent mixing, increases in case. When 
convection is very strong the convective growth rate becomes of order 0, and 

F r 1 

t££IE. ^ ^ _^ (57) 

Fturb a ' 

so that a small a implies a small role for turbulent mixing. 

Finally, we note that the results of our model differs from previous work in two interesting 
ways. Firstly, the intermediate stable branch is completely suppressed and this implies that DNe 
accretion disks do not come to rest, even briefly, midway through the ionization transition. This 
result would weaken the "stagnation effect" ( Mineshige 1988| ) which relied partly on the existence 



of this middle unstable branch to account for UV delay in the early rise of DN outbursts. One 
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would still have stagnation, however, due to the large specific heat at the temperatures between 
~ lO'' and 2 x 10^ K. Secondly, the value of T,max is reduced by about 20% and this is quantitatively 
interesting. However, we see from this figure that this correction is less important than arriving 
at a definitive value of ao in the scaling law. Choosing oq = 20, which is a conceivable, albeit 
extreme, choice, can almost double T,max while producing only a minor change in Timin- 

Fig. ^ shows the fractional heat fluxes as a function of height within a DN accretion disk. 
The top panel shows a model with turbulent transport and damped convection while the bottom 
one is a model with undamped convection. These two cases have the same total half-thickness 
log /i(cm) = 8.80 and lie on the middle branch of the S-curve. The undamped convection model 
is roughly at the upper end of the intermediate branch. Both models have log r(cm) = 10.5 
and use the scaling law for a. There is no significant change in F^ad between the two models, 
indicating that the total mixing flux is roughly the same, but in the case of damped convection 
this flux is almost entirely carried by Fturb- We note that formally the convectively unstable region 
(V_Rc > ^ ad) is not dramatically smaller in the damped case. Instead, while both models have 
broad, and roughly equal, convectively unstable zones (~ h), the shrunken convective zone shown 
in Fig. |2| in the damped case is due to the lack of solutions to the minimization problem for weak 
convection. This result for weak convection shows that the presence of an available source of free 
energy is not a guarantee of instability if the environment supplies strong damping. In Table |l] we 
see that where convection is strongest in the undamped convection model, it just barely survives 
in the damped model, even though the eddies are almost isotropic and roughly adiabatic (F' > 1). 
At this level, even their complete elimination would not change the S-curve significantly since 
the curve for Fau,damp already coincides with the convectionless model. The factor 1/(1 -|- ooF") 
describing the turbulent losses is smaller as F' increases. This means that convective and turbulent 
mixing becomes stronger (or weaker) than radiation simultaneously and this in turn explains why 
turbulent mixing is not important on high and low temperature stable branches of the S-curves 
where radiative transport dominates and the convectively unstable region is narrow. Although 
Fturb exists in the convectively stable region and carries energy toward the mid-plane, increasing 
the temperature gradient, Ft^rb behaves in other respects like Fconv with a shorter mixing length 
and is severely damped by radiative losses when radiative transport is important. Our calculations 
show that F" ~ 10^^ in the stable hot and cold states for the Fj-ad + Fturb curve in Fig. |l|, 
rendering mixing transport of any kind largely irrelevant. This is why all the S-curves coalesce in 
the hot and cold states. Fig. ^ also shows that for both models turbulent mixing dominates near 
the midplane where the disk is more opaque and where convection is weak due to small Qz- The 
numbers in parentheses in Table |l| (or in Table |^, and Table ||) are the corresponding results for 
the case F^ii nodamp for the same h. As expected, these numbers are larger than those for damped 
convection. The growth rate, 7, for undamped (MLT) convection in these tables is Fc/\/2 as 
defined in the preceding section. 

Table § shows the values of a along the S-curves in Fig. |l[ The outburst duration for DNe 
seems to imply that a ~ 0.2 in the hot state (cf. gmak 1984 

annizzo 19941 ). Evidently the 
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scaling law leads to a considerable overestimate (up to about 0.6) as we go up the hot state branch 
of the S-curve for DN disks. In addition, the use of the scaling law for a in the hot state in DN 
disks leads to a contradiction with the Bailey relation which says that the decay time associated 
with the decline from outburst in DNe increases with orbital period. Both of these points suggest 
that the scaling law should start to saturate as a rises above ~ 0.1. We note that the scaling law 
is motivated by the behavior of SXTs, which have no Bailey relation and much stronger evidence 
for an exponential decline from outburst. One simple compromise might be to let a rise to a 
maximum value of 0.2 and then abruptly go over to a constant. This would have no effect on 
the middle and lower branches of the S-curve, where turbulent mixing and the suppression of 
convection can have some effect. However, we note that a more physically motivated geometrical 
dependence for a involves a prolonged transition from a steep scaling to essentially constant 
behavior ( [Vishniac 1999| ). An alternative possibility is to choose a smaller value of ao, equivalent 
to asserting that the central black hole mass in SXT's is smaller. Table ^ also displays the values 
of a based on the formula 20(cs/r/ri)^/^ along the S-curve indicated by the dash-dotted line in 
Fig. |l], equivalent to asserting that the central mass in SXT systems is typically ~ 4Mq, which is 
a very low value. 

The S-curves for SXT disks are illustrated in Fig. ^ at r = lO^'^'^ cm with the central object 
of 10 Mq. The much larg er Yjjjiin and ^max along with the more pronounced hysteresis curve 
for the total flux {Frad + Fconv,MLT) reflect on stronger convection in SXTs than that in DNe at 
log r (cm) = 10.5. This results from the larger S for SXTs. The energy flux radiated from the 
surface is given by 

IVm'^ ~ i/SJ]^ ~ ac^SO. (58) 

Since S-curves for different systems operate at the similar T^ff, k (i.e. the condition for partially 
ionized hydrogen) and oc Tmid oc r^/^ oc S^/^, the above equation implies 



1 

. (59) 



Moreover, the scaling law for a gives 



a oc 



) (60) 



MJ 



which actually does not depend strongly on S. Together with equation (59) this gives 

S oc M4/23ri2/23^ (61) 

in accordance with the usual result that the column densities in the ionization transition region 



drop as a cooling wave propagates inward in a disk (for recent calculations see Cannizzo, Chen, 



and Livio 1995, Hameury et al. 19981). Similarly, we can obtain a crude estimate of the strength of 



convection on the unstable branch from equation (p^). Since k is roughly the same for disks with 
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similar ionization fractions and temperatures, and since we have already seen that the midplane 
temperature differences will be slight, we can rewrite equation (17) as 

i'conv,MLT/rrad^ ^ CX2.i2 ~2T746~' ^ ' 

Cs T I 

where we have ignored the change in vmlt/cs- This implies that convection is relatively more 
important for larger central masses and at smaller radii. In fact, SXT disks have larger column 
densities than DN disks for ionization transitions at the same r and consequently stronger mixing. 
Fig. ^ shows that the curves {F^ad + Fconv,MLT) and (-F^ad + Fturb) are more different than in 
Fig. ^. Conversely, the smaller difference between the curves Fau,nodamp and {F^ad + Fconv) also 
illustrates the increased importance of convection (and reduced importance of turbulent mixing). 
Adding to this effect is the smaller a which reduces the importance of turbulent mixing and 
damping. The end result is that the damped F^onv in SXTs is stronger than in DNe. Although 
the Faii^damp curve is still close to the curve with no convection, it consistently lies to the left and 
includes a very small intermediate branch. As before the value of ^max is reduced, in this case by 
about 30%, relative to the old standard model. Following Table ^ we note that the values of a 
used for these models never reach the unrealistically high values seen on the hot branch in DNe. 

Fig. ^ is corresponds to Fig. but for SXTs instead of DNe. In this case the total height for 
both models is log h (cm) = 8.33. The damped convection model in the top panel lies near the 
upper end of the short intermediate branch. As expected, convection is much more significant in 
SXTs in terms of both its intensity and its vertical range. This is also seen in Table |3| which shows 
that Ptot is of the order of unity in the confined convective region, suggesting that damping is 
weak for the strongest convective modes. In addition, the difference between that of damped and 
undamped convective speeds is much smaller than for the DNe cases shown in Table |l[ Table |^ and 
Table ^ (including Table || which we will discuss later in AGN cases) also show that for sufficiently 
strong convection the aspect ratio A usually increases with height since radiative damping (F' ^ 1) 
and turbulent drag {k'^i', see equation (0)) are more severe farther from the midplane. In fact, for 
other values of h we find that slightly radially distended convective eddies {A slightly larger than 
1) can actually be favored at high altitudes. The survival of strong convection in the middle of 
convectively unstable zones and the complete suppression of weak convection once the convective 
growth rate drops raises the question of convective overshooting. The error that results from 
ignoring overshooting should be small in our case for two reasons. First, the total mixing flux does 
not change if damping is added. Second, overshooting is important in cases where no other mixing 
mechanism exists. In accretion disks there is always turbulent mixing. Finally, we note that in our 
case convective overshooting still carries energy away from the mid-plane in the overshoot region 
since Vrc > '^ad- This is in contrast to the stellar case, where Fconv is negative in overshoot 
region, which is convectively stable, and this in turn increases the temperature gradient. 

Fig. ^ shows the S-curves for SXT disks for log r (cm) = 9 and 9.75. Although the typical 
S is smaller than at a radius of log r (cm) = 10.5, as expected from equation (|6l|), the larger 
value of rJS^ (cf. equation (^)) and the smaller value of a (cf. equation ( |60[ ) and Table ^) make 
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convection slightly stronger and turbulent mixing slightly weaker. Consequently, the two curves 
,mlt) and Frad,nodamp coalesce, while we see a larger difference between the curves 
for {Frad ~l~ Ff^onv ,mlt) and (Frad + Fturb)- However, when we compare to Q, we see that the 
strength of convection drops at smaller Teff and that this effect is stronger at smaller radii. The 
end result is that the final curve, Faii^damp follows the convective curves at high temperatures and 
switches abruptly to the turbulent mixing curve at low temperatures when the aspect ratio of the 
convective cells drops to the point where radiative losses and turbulent mixing turn off convection. 
This leads to a pronounced intermediate branch at small radii. For log (r (cm) = 9 the reduction 
in Timax is less than 20%. 

Fig. P displays the S-curves for the usual AGN parameters: Mbh = 10^ Mq and 
log r (cm) = 15.5. AGN are extreme cases relative to DNe because of their much smaller a, 
which implies weaker turbulent damping and larger Fconv/ Fturb, and much larger S, hence larger 
Fconv/Frad- Consequently the curve Faii^nodamp hes close to the curve {Frad + Fconv,MLT) and both 
are far from the curve (F^ad + Fturb)- There is almost no difference between the curve Faii^damp and 
the curve Faii nodamp- of this is as expected for a disk where convection is strong and turbulent 
mixing is weak. The weakness of turbulent mixing is also revealed in the relatively large difference 
between T^id and T^jj for the Fj-ad + Fturb curves compared to the convection S-curves; that is, 
the temperature difference is increased in the middle branch as a result of inefficient turbulent 
mixing. We also ran the code for {Frad + Fconv,MLT) with a reduced mixing-length A = 0.316Ap 
with the usual scaling law for a. The S-curve in this case, depicted by the line composed crosses 
in Fig. y, is more or less similar to the CCHP model, which attempts to account for shear through 



a general reduction in the mixing length, at least for 0.01 ;^ a ^ 1 (Cannizzo 1992). We see 
that this approach exaggerates the suppression of convection for AGN. In our model convection 
remains relatively efficient for AGN despite shearing and turbulent damping. 

The relative importance of the vertical heat transport mechanisms for AGN is illustrated in 
Fig. 1^. We choose models with log h (cm) = 12.58 lying on the unstable branch. We see that the 
overall structure of the disk is affected only slightly by including turbulent damping of convective 
cells. In both cases turbulent mixing dominates near the midplane and decreases rapidly in 
importance away from it. Radiative transport plays a major role only in the region just below 
the photosphere. The vertical structure of the damped convection zone is also shown in Table |5[ 
We see larger convective efficiency factors associated with radiative losses P', small turbulent 
drag k'^u, and very small turbulent loss factors 1/(1 -|- aoT") compared to those in Table |l| and 
Table ^ As a result, Ptot is close to unity and convection is well approximated, at least in terms 
of geometry, by isotropic convection without turbulent mixing. As in the SXT cases, we also 
find that A is occasionally bigger than 1 in the most robust convective modes. The larger role 
for turbulent losses 1/(1 -|- aoF") for AGNs is really a consequence of very small radiative losses, 
and does not reflect an overwhelming suppression of convection by turbulent mixing. We can see 
this by considering the total efficiency factors obtained by dividing the F' without parentheses by 
the corresponding 1 -|- aoT" . In fact. Tables |l[ |^, and [s] show that the total convective efficiency 
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r'/(l + oor") for these three systems are of the same order. This is because the smah values of a 
and huge values of T, for AGNs assure that radial turbulent mixing and radiative losses are almost 
negligible for them while the turbulent flux is extremely weak. 

The very large values of S implied by the scaling law for a applied to AGN disks (see Table |6|) 
leads to the question of disk self-gravity. If the self-gravity of the disk is as important as the 
gravity imposed by the central object, then local gravitational instabilities may be important and 
our description of AGN disks is not self-consistent. In Fig. |^ we see the ratio between self-gravity, 
Qseif = 27rGS, and Qcentrai = for complete disk models of AGN with turbulent damping. We 
see that the two become comparable when the AGN disk Te/j approaches 10^ K. This should raise 
a and turn the curves sharply to the left. 



Fig. y and Fig. 10 illustrate the overall behavior of convective flux along the unstable branches 
of different accretion disk systems. The ranges of effective temperature displayed in these figures 
are associated with those for the thick solid lines (the curves Fau^amp) in Fig- ^ Fig. ^, and Fig. |^. 
Maximum values of fractional convective fluxes (both MLT and damped) and the vertical extent 
of the convective zones, Iconv (scaled with half-thickness /i), are depicted for damped convection as 
a function of T(,ff in Fig. ^. MLT convection (dashed lines) becomes stronger and stronger as one 
goes from DNe (top panel), to SXTs (middle panel), to AGN (bottom panel) and the difference 
between damped and MLT convection diminishes over the same range. The importance of damped 
convective flux can also be measured by Iconv (dotted lines), where Icxmv is the vertical range for 
which solutions to the minimization problem exist (for which convection can actually occur). 
On average, Iconv increases from 0.1 for DNe to ^ 0.8 for AGN. Fig. |1C| displays the maximal 
and minimal values of \/~A at a given Tcjj associated with the S-curves for Faii^damp (denoted 
by thick solid lines in figures |l|, ^, and ^ for different systems. As mentioned above, A tends 
to increase with height and so a larger separation between maximal and minimal A goes with a 
larger vertical range for damped convection. Values of y/A as small as ~ 0.3 can exist for AGN 
but not for DNe in accord with the result that convection is strong in AGN. In fact, Table |l| and 
Table |5| show that the convective growth rate 7 k, 0.1 for DN at ^/A « 0.8 is comparable to that 
for AGN at \/~A ~ 0.3. No solutions exist on the left side of dotted lines for weak convection as a 
result of the severe damping of narrower convective cells. On the other hand, 'fatter' {y/~A > 1) 
convective eddies can exist for strong convection, since these solutions reduce radiative losses and 
radial mixing and extract only a small penalty in the form of increased vulnerability to secondary 
instabilities. 

Metals are the main supply of free electrons in the cold, neutral accretion disk state. In our 
code the ratio of electron gas pressure to the total gas Pe/Pg is ~ 10~^ around Umax, indicating 
that the disks are only weakly ionized. We calculated the microscopic resistivities due to the 
electron-neutral collisions using the equation rf = 230(n„/ne)T^/^ cm^ s^^ (Blaes and Balbus 1994) 



and used the standard formula for electron-ion collisions, r/ = 2.65 x lO^^T ^/^cm^/s. We found 



- 23 - 



magnetic Reynolds numbers f\ Rcm = ^pCs/f] in quiescent DN and SXT disks are less than 10^, 
in accord with previous work ( |Gammie and Menou 1998| ), although this neglects the possibility 



of ionizing particles from the disk coronae. Fig. 11 illustrates how the S-curves for disks with 



solar metal abundances differ from ones without any metals around S^ax for different systems. In 
order to distinguish the effects of MLT convection from turbulent mixing, we use the curves for 
{Frad + Fconv,MLT) and (Frad + Fturb)- For curves with MLT convection, the discrepancy between 
metal and no metal models near cases around Smaa; increases as one goes from DNe, to SXTs, 
to AGNs. This is because convective mixing becomes stronger in AGNs. On the other hand, 
the opposite trend happens to the curves with turbulent mixing which follows from the fact that 
turbulent mixing becomes weaker in AGN. 

Finally, in Fig. (||) we test the consequences of using a constant value of a instead of the 
scaling law. This is motivated by Gammie and Menou's suggestion ( |Gammie and Menou 1998 , 



see also Armitage, Livio, and Pringle 1996) that the timing of DN outbursts can be explained by 
the relatively high resistivity of quiescent DN accretion disks. If the BH instability shuts off below 
ReM ~ lO"^, then the low values of a inferred from the intervals between outbursts reflect simply 
the operation of some different angular momentum transport mechanism in quiescent disks, and 



not a general scaling law. Fig. (12) shows the S-curve down to relatively low temperatures, and 
the associated range of Rcm- The important feature here is that a cutoff of Rcm ~ 10^ or less 
results in S-curves with a very small range in column density. While S^ax may still be very large 
if a drops sharply at still smaller temperatures, the propagation of a cooling front in an accretion 
disk depends to a large extent on the nature of the S -curve near S^^n and detailed models of 
cooling fronts indicate that the cooling sections of the disk trace a path just below the unstable 
equilibrium curve ( [Menou, Hameury, &: Stehle 1999| ). A cooling wave with a nearly vertical 



drop, as seen here, cannot be described by the rapid cooling wave theory (Vishniac and Wheeler 



1996 , Vishniac 1997| ) but should instead propagate at a fraction of the viscous accretion speed, just 



as expected for an S-curve with a small value of 'Smax/'^min ( [Vishniac 19981 ). In fact, models with 
constant a not only show slow cooling front speeds but also frequent reversals ( ^mak 1984 , Menou, 



Hameury, &: Stehle 1999| ). Such models lack a sharp drop in a at very low temperatures and may 



not represent realistic calculations of the results of Gammie and Menou's model. However, it is 
not clear that even the existence of a second unstable branch at very low temperatures would 
affect the cycling of a DN disk. The S-curve shown here is sufficient to produce a complete, albeit 
weak, thermal cycle without ever moving to lower temperatures. We caution however that this 
conclusions depends on using a particular value of Rcm for a cutoff. If the cutoff value of Reu lies 
in the range 10^ to 10^ then a robust thermal limit cycle is still be possible for a constant a model. 



^The more formal definition of Rcm for the isotropic turbulence driven by BH instability is i^/rj. 
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7. Summary and Conclusions 

We have proposed a new model for vertical energy transport in optically thick, Keplerian 
accretion disks which takes into account the effects of turbulent mixing due to the BH instability. 
This includes both an evaluation of the turbulent flux Ff^rij driven by approximately isotropic 
turbulent cells with an effective vertical thermal conductivity ~ u, and the increased damping of 
the sheared convective cells. The latter effect includes not only radiative losses, which have been 
broadly discussed in the literature, but also turbulent mixing which enhances heat transport and 
retards convective motion. Simple physical arguments, and the linear dispersion relation, indicate 
that weak (Tc < Jl) convective eddies are radially thin in the shearing flow of an accretion disks. 
This enhances the effects of turbulent drag, radial mixing, and radiative losses. We account for 
these effects through an extension of linear perturbation theory using a one mode model. We 
choose the linear mode most resistant to secondary instabilities and to turbulent damping and 
by calibrate the effects of this mode by requiring agreement with MLT if turbulent viscosity and 
shearing are removed. Following this procedure, our model for convection becomes a minimization 
problem for the function Ptot subject to the linear dispersion relation involving several damping 
mechanisms. 

Our use of quasilinear theory, calculating quadratic transport quantities and typical eddy 
shapes using linear modes, is reasonable, but not rigorously justifiable. Generally this procedure 
works only when energy flows through a system in ways that already appear in a linear analysis. 
The appearance of strong finite amplitude instabilities would destroy the rationale behind this 



approach. Fortunately thin accretion disks do not appear to possess such instabilities (Hawley, 
Balbus and Winters 1999). One of the other major components of our approach, the suppression 
of one family of unstable modes (convection, in this case) by another is unusual, but not 
unprecedented. In laboratory plasmas we have the example of radial convection being disrupted 
by the spontaneous appearance of shear (cf. Burrell 1997). In the context of accretion disks 
the suppression of the Parker instability by BH induced turbulence ( ^tone et. al. 1996| ), and 
previously predicted in Vishniac and Diamond (1992), is very closely analogous to suppression of 
convection. The most problematic aspect of our treatment is our method for choosing a 'typical' 
linear mode (cf. equation (|3^). As far as we know this method has never been proposed before 
and there are no definitive tests of its usefulness. However, it does capture the correct limits, i.e. 
settling on the remaining unstable mode in the case of marginal stability and choosing the modes 
least susceptible to secondary Kelvin-Helmholtz instabilities when damping is unimportant. We 
conclude that it is at least qualitatively correct. 

The thermal equilibrium S-curves calculated in our model differ from previous models in 
several ways. Our models tend to reproduce previous work along the hot stable branch of the 
S-curve, and give the same, or nearly the same Smm- However, the unstable branch tends to lie 
somewhat below, and to the left, of previous work in the T, — T plane. In addition, the short, 
intermediate stable branch seen in previous work often disappears, and is usually more nearly 
vertical when it does appear. The value of Smax is smaller in our models, though only by 20-30%. 
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Crudely speaking, these effects are stronger for thicker disks, that is, for DNe and SXTs, rather 
than AGN. However, in some respects these effects are more dramatic for SXTs than for DNe, 
largely because the different transport mechanisms give rise to more divergent S-curves in the 
former. AGN disks are still reasonably well described by standard MLT theory without turbulent 
mixing. There are two reasons for the way turbulent effects scale from DNe to AGN. First, in an 
average sense, convection becomes stronger as we go to thinner disks around more massive objects. 
Fconv / Frad rJS^. This leads to faster convective growth rates, which means that convective cells 
become more resistant to disruption by any reasonable level of local turbulence. Second, we have 
adopted a scaling law which reinforces this trend by lowering the value of a, and consequently 
the strength of turbulent dissipation, in thinner disks. Since a decreases and S increases as one 
progresses from DNe, to SXTs, to AGNs, we find an almost complete suppression of convection in 
DNe. Models with no convection are reasonably good approximations to our models with damped 
convection. 

The specific calculations presented here rely to large extent on rather uncertain approximations 
in extending MLT to strongly sheared environments and in estimating the level turbulent dissipation 
in disks. While the procedure followed here is reasonable enough, it would clearly be better to 
have means of calibrating our theory for anisotropic convection. Nevertheless, the qualitative 
nature of our results are not sensitive to the specific choices we have made. On the whole, the 
results are less sensitive to the exact level at which weak convection is suppressed than to the 
qualitative point that weak convection is vulnerable to turbulent mixing in a sheared environment. 

The work presented here is based on the assumption that the scaling law equation (^) for 
a is a reasonable approximation. This tends to make our corrections to MLT less important for 
thin disks, especially AGN disks. In fact, the status of this scaling law is somewhat uncertain, 
for reasons given above. However, we note that a constant q, down to some limiting value of 
the magnetic Reynolds number, seems unable to reproduce the basic phenomenology of DNe 
outbursts. At the same time a very low a seems required from observations of quiescent DNe (for 
example Wood et al. 1986, 1989| ). Unless the cooling wave theory is discarded altogther it seems 



necessary to postulate some sort of geometric dependence for a. On the other hand, the Bailey 
relation for DNe suggests that the scaling law saturates for values of a typical of DNe in outburst. 
Finally, we note that the role of irradiation ( [King fc Ritter 1998 ) and possible advection flows 



( Narayan, McClintock, and Yi 1996 ) in SXTs remain unclear, and prevent a clear understanding 



of the form of a required to explain observations of thin disks. 
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Table 1. The Structure of Damped Convection in DN disks: log /i (cm) = 8.80, 
log r (cm) = 10.5, a = 8.89 x 10^^ 



z 


r' 


7/n 


Va 


Ptot 


v/v 






1/(1 + aor") 


0.67 


2.01(5.01) 


0.07(0.51) 


0.81 


117 


1.10(2.07) 


0.03(0.21) 


0.76 


0.06 


0.70 


1.72(3.29) 


0.12(0.56) 


0.88 


45 


1.16(1.98) 


0.05(0.22) 


0.77 


0.11 


0.72 


1.02(1.57) 


0.16(0.61) 


0.95 


26 


1.18(1.70) 


0.06(0.23) 


0.81 


0.21 


0.74 


0.12(1.57) 


0.06(0.58) 


0.91 


249 


1.04(1.44) 


0.02(0.21) 


0.99 


0.42 



Table 2. Values of a in the upper, middle, and lower branches of S-curves for DN disks in Fig. |T|. 



QO 


branch 


a 


50 


upper 


^ 0.25 - 0.61 




middle 


^ 0.02 - 0.25 




lower 


< 0.02 


20 


upper 


0.11 - 0.26 




middle 


^ 0.007- 0.11 




lower 


< 0.007 



Table 3. The Structure of Damped Convection in LMXB disks: log hem = 8.33, log rem = 10.5, 

a = 1.55 X 10~2. 



z 


r' 


7/n 


Va 


Ptot 


V/V 






1/(1 + aor") 


0.47 


34.5(83.5) 


0.07(0.21) 


0.55 


23 


1.25(1.51) 


0.03(0.09) 


0.30 


0.01 


0.66 


14.4(15.7) 


0.18(0.36) 


0.64 


3.2 


1.49(1.50) 


0.07(0.14) 


0.21 


0.09 


0.72 


6.63(6.95) 


0.25(0.44) 


0.72 


2.4 


1.53(1.52) 


0.09(0.15) 


0.22 


0.20 


0.79 


1.11(1.35) 


0.33(0.52) 


0.81 


2.0 


1.35(1.36) 


0.10(0.15) 


0.28 


0.61 
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Table 4. Values of a in the upper, middle, and lower branches of S-curves for SXT disks in 

Fig. I and Fig. H 



ao 


branch 


a at log r (cm) — 10.5 


a at log r (cm) — 9 


50 


upper 


^0.06-0.15 


Z 0-004 




middle 


^ 0.003 - 0.06 


« 0.0002 - 0.004 




lower 


^ 0.003 


^ 0.0002 



Table 5. The Structure of Damped Convection in AGN disks: log h (cm) = 12.58, 

log r (cm) = 15.5, a = 8.32 x 10"''. 



2 


r' 


7/^7 


Va 


Ptot 


V/V 


Vz/Cs 




1/(1 + aor") 


0.28 


7.73x10^(7.40x10®) 


0.05(0.11) 


0.28 


10.8 


1.06(1.03) 


0.01(0.03) 


9.92x10"^ 


1.70x10"'' 


0.46 


4.91x10^(2.54x10'') 


0.11(0.16) 


0.35 


4.6 


1.07(1.05) 


0.03(0.04) 


9.45x10"^ 


5.61x10"'' 


0.66 


9.33xl0'*(2.46xl0^) 


0.28(0.41) 


0.56 


2.1 


1.15(1.12) 


0.05(0.06) 


8.81x10"^ 


5.94x10"^ 


0.85 


1.66x10^(3.05x10=*) 


0.93(1.12) 


0.90 


1.13 


1.60(1.56) 


0.11(0.12) 


9.66x10"^ 


7.90x10"=* 


0.93 


1.63(2.84) 


2.99(3.52) 


0.99 


1.02 


1.80(2.25) 


0.24(0.27) 


1.68x10-1 


9.36x10"^ 



Table 6. Values of a in the upper, middle, and lower branches of S-curves for AGN disks in 

Fig. |. 



50 80000 - 210000 ^ 0.004 - 0.02 

37000 - 80000 ^ 0.0018 - 0.004 

14000 - 37000 ^ 0.0006 - 0.0018 

5000 - 14000 ^ 0.0002 - 0.0006 
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200 400 600 800 

i;(g/ cm^) 

Fig. 1. — S-curves of DN disks {Mwd = ^Mq, log rem = 10.5) for different vertical energy 

transport models: Faii^amp (thicker solid line), Faii,nodamp (dotted line), Frad + Fconv,MLT (dashed 
line), and F^ad + Fturb (solid line). While these curves are based on the scaling law a = 
50(cs/r/r2)^/^, the dash-dotted line denotes the S-curve for Fj-ad + -^turb based on the scaling law 
a = 20(c,/r/J7)3/2. 
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Fig. 2. — Fractional energy fluxes for DN disks (M^i/d = l-^O) log rem = 10.5) along the vertical 
direction z for the same half-thickness log /i cm = 8.80 in the cases of damped convection (top 
panel) and MLT convection (bottom panel). Frad/Ftotai, Fconv/ Ftotal, and Fturb/Ftotal are denoted 
by a solid line, a dashed line, and a dotted line respectively. 
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Fig. 3. — S-curves of SXT disks {M\yd = 10 Mq, log rem = 10.5) for different vertical energy 
transport models: Faii,damp (thicker solid line), Faii,nodamp (dotted line), Frad + Fconv,MLT (dashed 
line), and Frad + Fturb (solid line). All curves are based on the scaling law a = 50(cs/r/r2)^/^. 
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Fig. 4. — Fractional energy fluxes of SXT disks = 10 M©, logr cm = 10.5) along the vertical 

direction z for the same half-thickness log /i cm = 8.33 in the cases of damped convection (top 
panel) and MLT convection (bottom panel). Frad/FtotaU Fconv/Ftotah and Fturb/Ftotal are denoted 
by a solid line, a dashed line, and a dotted line respectively. 
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Fig. 5. — S-curves of SXT disks around a 10 central black hole at log rem = 9 and 9.75 for 
different vertical energy transport models: Faii^damp (thicker solid line), Faii^nodarnp (dotted line), 
Frad + Fconv,MLT (dashed line), and Fpad + Fturb (solid line). All curves are based on the scaling 
law a = 50(cs/r/ri)^/^.. We also show the aspect ratio of convective cells, A^^"^ and the strength 
of convection Fc/ri in the middle of the convective region as a function of Tgff. The values of A^/'^ 
and Fc/J^ are the maximum and minimum values for a given T^ff- 
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i;(g/ cm^) 

Fig. 6. — S-curves of AGN disks {Mwd = 10^ Mq, log rem = 15.5) for different vertical energy 

transports: Faii,damp (thicker solid line), Faii^nodamp (dotted line), Frad + Fconv,MLT (dashed line), 
and Frad + -^twrb (solid line). The line composed of crosses denotes the S-curve for the flux model 
Frad + Fconv,MLT with a shortcr mixing-length log(A/Ap) = —0.5. All curves are based on the 
scaling law a = 50(cs/r/r2)^/^. 
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Fig. 7. Fractional energy fluxes of AGN disks [My/D = 10^ M©, log rem = 15.5) along the 
vertical direction z for the same half-thickness log ft. cm = 12.58 in the cases of damped convection 
(top panel) and MLT convection (bottom panel). Frad/Ftotal, Fconv / Ftotal , and Fturb/Ftotal are 
denoted by solid line, dashed line, and dotted line respectively. 
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Fig. 8. — Effective temperature T^ff vs. the ratio of self-gravity to central object gravity 
9self / Qcentral for AGN models whose total vertical energy flux is Faii^amp- 
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Fig. 9. — Maximum values of convective fluxes and vertical ranges of damped convection are 
plotted with respect to effective temperature Tgff for different systems. Tlie dashed lines denote 
Fconv,max/Pall,nodamp, where Fconv,max is the maximum MLT convective flux at a given Teff. The 
solid lines denote Fconv,max/Fall,damp, where Fconv,max is the maximum damped convective flux at 
a given T^ff- The vertical ranges (scaled with h) of damped convection are indicated by the dotted 
lines. 
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Fig. 10. — Effective temperature vs. the ratio of radial to vertical length ^fA of damped convection 
for different systems. The solid and dotted lines represent the maximal and minimal values of \fA 
respectively at a given Tg//. No solutions exist on the left side of dotted lines for weak convection 
(i.e. smaller \fA) as a result of strong damping. 
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Fig. 11. — Influence of metal abundance on the S-curves around Umax for MLT convection and 
turbulent mixing in DN (top panel), SXT (middle panel), and AGN (bottom panel) disks. Solid 
lines indicate the curves Ffot = F^ad + Fconv,MLT with metals. Dashed-dotted lines indicate the 
curves Ftot = Fj-ad + FcouvMLT without metals. Dashed lines denote the curves Ftot = Frad "I" -^turb 
with metals. And dotted lines denote the curves Ftot = Frad + Fturb without metals. 
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Fig. 12. — S-curves of DN disks {Mwd = Mq, log r(cm) = 9) associated with different energy 
transports for a = 0.1 around ^min- Faii^amp (thicker sohd hne), Faii^nodamp (dotted hne), 
Frad + Fconv,MLT (dashed line), and F^^d + Fturb (solid line). 



